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ABSTRACT 



Optimal age replacement policies are designed to cut down system failures and 
minimize maintenance cost. By scheduling planned replacements, a system is replaced 
at age 0* or at failure, whichever comes first, and the cost of replacement before failure 
(planned) is less than the cost after failure (uiq>lanned). In this thesis, the distribution of 
lifetimes is a known, increasing failure rate phase type distribution. To find the optimal 
age of replacement, the parameters of the underlying phase type distribution must be 
estimated. 

An optimal age sequential estimation procedure is developed. In particular, the 
phase type distributions parameters are estimated using a matching moments nonlinear 
programming approach. Since there are many parameters associated with phase type 
distributions and the distributions include matrix exponential terms, the parameters are 
in general difficult to estimate. A specific case where the phase type distribution has 
initial probability vector a=(l,0,0) is studied for different sample sizes and compared 
with a similar nonparametric procedure. 
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I. INTRODUCTION 



Normally a system or component is replaced when it fails and Cj, the cost of 
replacement before failure (planned replacement), is often less than C 2 , the cost after 
failure (ui^lanned replacement). The problem of determining when to rq)air and when 
to replace failing systems is the concern of management resources. Inefficient 
management due to the use of nonoptimal maintenance policies can lead to significant 
system maintenance cost. In general, optimal maintenance policies are designed to cut 
down the number of system failures and minimize maintenance costs. In this thesis, we 
study the problem of estimating a type of optimal maintenance policy (the optimal age 
replacement policy) when the underlying system life distribution is phase-type. In 
particular, we consider an adaptive estimation scheme where the estimated optimal 
replacement age is updated each time the system is replaced. 

Maintenance is defined to be all activities taken to keep the system in serviceable 
condition or to bring it back to serviceability. There are two types of maintenance, 
corrective and preventive maintenance [Ref. 5: pp. 1-8]. Corrective maintenance is used 
after a failure. This does not necessarily mean that such action has not been foreseen. 
Preventive maintenance aims to reduce the probability of failure and includes: 

-Planned (scheduled) maintenance in which specified components are replaced 
(usually at regular intervals). The maintenance time is usually based on component 
lifetime distributions. 
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-Unplanned (condition-based) maintenance in which the decision to replace or not 
to replace is made according to the outcome of a diagnostic study. 

The policies that are designed to reduce the number or the probability of system 
failure and maintenance costs are called maintenance policies [Ref. 21: pp. 1-3]. We 
concentrate on age replacement policies where a system is replaced at age T or failure, 
whichever comes first. For this problem to make sense, Q must be less than Cj and the 
system must age with time i.e. the failure rate of the system must be increasing with age. 
It is not wise to replace the equipment before failure when the system has a constant or 
decreasing failure rate such as when the system failures occur according to an exponential 
distribution. In this case, replacement will not reduce the probability of system failure 
occurring in the next instance of time. When using an age replacement policy the 
question always asked is, "At what age should the equipment be replaced". If the 
replacement occurs too frequently the maintenance costs will be high. If replacement is 
too infrequent, the system will fail more often than necessary and again, the maintenance 
cost will be too high. Thus it is desirable to find an "optimal" replacement age. Here, 
the optimal age is defined as the age which yields the minimum long run expected 
maintenance costs per unit time. 

We will assume that the maintenance action returns the system to the good-as-new 
condition, thus the same services are provided as before replacement. To accomplish this 
we assume that the systems used for replacement have lifetimes that are independent and 
identically distributed according to a phase-type distribution. The class of phase type 
distributions is large and includes Exponential and Gamma distributions along with 



2 



convolutions and mixtures of these distributions. This fact, along with the fact that phase 
type distributions have a physical interpretation make them particularly well suited for 
modeling system lifetimes. 

The optimal replacement age depends on the underlying phase type distribution. 
This is in general unknown and must be estimated. Estimation for phase type 
distributions based on iid lifetimes has been studied by Neuts and Meier (1980). 
Estimation of the optimal age of replacement for phase type distributions is new. 
However both parametric and nonparametiic estimation based on iid observations have 
been considered by Arunkumar (1972), Ingram and Schaeffer (1976), Bergman (1977), 
and Barlow (1978). Here we use a sequential approach for estimating the optimal 
replacement age. At each replacement the estimate is updated and the new system is 
subject to an age replacement policy based on the optimal age estimated so far. This type 
of sequential approach has been studied parametrically by Oclay (1990) and 
nonparametricly by Bather (1977), Frees and Ruppert (1985), Aras and Whitaker (1991), 
Wu (1990) and in a decision theoretic framework by Glazebrook, Bailey and Whitaker 
(1991). 

Phase type distributions are described in Chapter n. The sequential estimation 
procedure to estimate the optimal age of replacement is given in Chapter HI, and in 
Chapter IV we present results comparing the sequential estimation procedure assuming 
an underlying phase type distribution with a nonparametric procedure. Conclusion and 
recommendation are given in Chapter V. 
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n. PHASE TYPE DISTRIBUTION 



A. GENERAL 

A phase type distribution is derined as the distribution of the time until absoiption 
in a finite state Markov process. To determine the distribution of the absoiption time, 
consider an m+1 state, continuous time Markov chain whose infinitesimal generator Q 
has the form 

^ It T 
||0 0 

where T is a nonsingular m X m matrix with negative diagonal elements and nonnegative 
off-diagonal elements, T° is vector of length m with nonnegative elements and 
0=(0,0,...,0). Moreover, 

Te -I- T° = 0' 
where e = (1,1,...,!)'. 

By the construction of Q the state m+1 is absorbing and all other states are 
transient. The necessary and sufficient condition for this is that the inverse R* exists 
[Ref. 15: pp. 446]. Let the initial probability vector of the Markov chain be given by ( 
o!,n+i ) with ae + = 1 and a being the m dimensional initial probability vector 

of transient states such that 0 < (ve 1 . Let X be the time until absorption into the 
(m+l)“ state. The probability distribution F(x) of the time until absorption in the state 
m+1 corresponding to the initial probability vector (a, ) is given by 
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F(x) = 1 - aexp(Tx)e , for x ^ 0 , 



( 2 . 1 ) 



where exp(A) is the matrix exponential of the matrix A defined in Ref. 8. 

The probability distribution F(x) on [0, oo) is said to be of phase type (PH- 
distribution) and the pair (a, T) is called a representation of F(x). If a„+i > 0 the 
distribution F(x) has a jump of height at x = 0. 

On (0, oo), the probability mass is distributed according to a density given by 



B. COST FUNCTION 

Let X, , Xj, . . . be a sequence of independent and identically distributed (iid) positive 
random variables with distribution function F. The sequence Xj, X 2 ,... represents the 
sequence of system lifetimes that would be observable if the system were replaced at 
failure. Let C,, CI 2 (Cj > Cj) and <l> be the respective cost of planned replacement; 
unplanned replacement and the replacement age. The long run expected cost per unit time 
can be verified [Ref. 1: pp. 87] to be 



where F(u) = 1 - F(u) is the survival function. A sufficient condition that guarantees a 



unique and finite <j>* that minimizes C(^), is that F have strictly increasing failure rate. 



f(x) = aexp(Tx)T° , for x > 0 , 



( 2 . 2 ) 



The noncentral moments n, = E[X‘] of F(x) are aU finite and given by 
/i; = (-l)4!(aT‘e) , for i 2: 1 [Ref. 15: pp. 446] . 



(2.3) 




( 2 . 4 ) 
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When the distribution F is phase type with representation (a, T), the long run expected 
cost function is 



C(<|)) 



Cj [1 -aexp (a\J>) e] + [aexp (2\|>) e] 




aexp ( Tu) edu 



(2.5) 



C 2 - aexp ( Jj>) e[C; -Cil 
[exp (2\|>) - J] e 



where I is an identity matrix. 

The phase type distribution does not necessarily have increasing failure rate, so the 
cost function C(<^) does not necessarily have a unique finite minimum. In the m =3 case 
even when the initial probability of the transient states is fixed the minimum can be 
unique and finite or can occur at infinity as shown in Figure 1, where a = (1, 0, 0) and 
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Figure 1 The long run expected cost per unit time curves of PH distribution with 
the same a and different T = T, or Tj 
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Figure 2 The long nin expected cost per unit time curves of PH distribution with 
the same T and different a 
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m. THE SEQUENTIAL ESTIMATION PROCEDURE 

In theory, any distribution of a nonnegative random variable can be approximated 
arbitrarily closely by a PH distribution [Ref. 11; pp. 1-2]. This means that the PH 
distributions are dense in where .^is the set of distributions with support on [0,<»). 
The paper of Johnson and Taaffe (1988) shows that the finite mixtures of Erlang 
distributions and PH distributions with a Coxian representation are both dense in Since 
both families are subsets of the family of PH distributions, this implies that PH 
distribution are dense in .^[Ref. 10: pp. 1-8]. Thus, the PH distribution can be used to 
approximate any unknown distribution of lifetimes (lifetime is always greater than or 
equal 0). Another feature that makes PH distributions a desirable choice to model system 
lifetimes is that they may have a physical interpretation. The absorbing state represents 
system failure and the transient states may represent different levels of a system’s ability 
to function. 

A. PARAMETER ESTIMATION 

There are many ways to estimate the parameters of a distribution including the 
method of maximum likelihood (MLE), and the method of moments. The PH 
distributions have special properties which make estimation difficult. First, the number 
of parameters is not fixed, it varies with the number of transient states, e.g. a PH 
distribution with 2 transient states has 6 parameters, a PH distribution with 3 transient 
states has 12 parameters etc. Also the probability distribution and density function include 
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the exponential of a matrix. This makes it hard to work with the likelihood function. In 
this thesis, the moment matching method is used to estimate the parameters of a phase 
type distribution but the MLE method and its associated problems are also discussed . 
In addition, the number of transient states, m, is assumed to be known. 

1 . The Maximum Likelihood Method 

Let X,, Xj,... represent the sequence of system lifetimes, where Xi, Xj,... 
are iid PH distributions with representation (a, T). After N observations the data 
available for estimation will be (Zj Si), (22,62), ...,(^,6 n) where Zi=min Xj is 

the i* lifetime and $\.t is the optimal age replacement estimated after i-1 observations and 
5 i=l if Xi^$*i_i and 6 j =0 otherwise. Assuming system lifetimes have a PH distribution 
the likelihood for the replacement ages Z,,...,Zn and types of replacement 6j, 62,...,6jj 
is 

jr 

L(a,T) = JJ [aexp(rz^)!T*]*^taexp(rZj)e]^‘*^ 

i=l 

N 

= JJ [-aexp(rz^)ire]*Maexp(rZi)e]^'*^ . 

1=1 

It is too difficult to maximize this likelihood directly by differentiating L(a,T) and 
solving for the MLE’s. There are several reasons for this. The likelihood and its 
derivatives include the exponential of a matrix form that cannot be written in closed form 
and the parameters are subject to numerous constraints. In addition, there are many 
likelihood equations due to the number of parameters (3 transient states will have up to 
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12 equations). An alternate approach is to And approximate MLE’s using nonlinear 
programming algorithms i.e.: 

If 

MAZL(a,T) = JJ [-aexp(rZi)re]*^[aexp(rz_£)e]^'*^ 

i=i 

S.T. a ^0 
ae = I 

tji < 0 for i = l,2,...,m 

tij ^0 for i j 

^ -tii fori = 1 , 2 ,..., m 

det(T) 0 

where m = the number of transient states of the PH distribution. Even when m is 
assumed to be known, there is still the problem of approximating the exponential of a 
matrix and the optimization software to take care of this problem is not available. Thus, 
it is very difficult to use this approach to estimate the PH distribution parameters. 

2. The Moment Matching Method. 

The PH distribution is a complicated distribution, there are many parameters 
and there are difficulties with computing the exponential of a matrix. One property, the 
existence of T*, implies that all noncentral moments of a PH distribution are finite. The 
k“* moment is given by 

= (-l)^!aT‘e. (3.1) 

Moment matching is a common method for estimating. Using moment matching by- 
passes the problem of evaluating the exponential of a matrix [Ref. 12: pp.3]. 
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Let m be the number of transient states of the PH distribution. There are m(m+l) 
parameters that need to be estimated, thus the m(m+ 1) equations from the moments that 
need to be solved are 

p,i = -aT’e 
p ,2 — 2!ofT"*e 



Mc(c+1) = (-ir^“^»K+m)!a'T‘e, 
where p,^_ is the k* sample moment. 

This method is still inappropriate due to the large number of equations and 
the fact that when the dimension of matrix T is big and the moments are of high order, 
substantial amount of error is introduced when trying to solve these equations. Johnson 
and Taaffe have shown that the moment matching method and nonlinear programming 
approaches can be used to estimate the parameters of a PH distribution. They match only 
the second and third standardized moments [Ref. 11: pp. 3-11]. 

Let Hi be the i* central moment. The second standardized moment, the coefficient 
of variation (C), defined as 

^ _ standard deviation 
mean 

can be written as 
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c 






(3.2) 



The third standardized moment, the coefficient of skewness, (7) is defined as 



Y 



^^3 



(3.3) 



And the t* standardized moment is for t = 3,4,.,, . 

The reason for matching standardized moments is that the standardized moments do not 
change with scale changes. So, the shape of a PH distribution with representation (a,T) 
will not change if it is multiplied by a nonnegative number. 

The nonlinear programming formulation for matching C and 7 to a 
continuous PH distribution is given by 

MIN. Wi(C-C(0))^ + W2 (7-7(0))^ 

S.T. 

oc ^ 0 
ae = \ 

tii < 0 for i = 1,2, ...,m 

tjj ^ 0 for i j 

fori = 1,2,.., ,m 

det(T) jii 0 

where 

C(0) = the second standardized moment of the current solution. 

7(0) = the third standardized moment of the current solution. 
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m = the number of transient states of the PH distribution. 



tjj = the elements of row i and column j of nonsingular matrix T. 

Wi,W 2 are positive weights chosen to guide the search and w, S: 3wj. 

The moments are matched when the objective function value is zero. The solutions which 
match the target moments are called "moment-matching solutions" [Ref. 11: pp. 5]. 

Even with fixed dimension, this procedure may have many solutions and these 
solutions are unpredictable. There are some practical points that need to be made 
associated with the properties of PH distribution and the problem of using nonlinear 
programming to get the moment matching solutions: 

- Different combinations of initial solutions and algorithm parameters may lead to 
different moment-matching solutions. Also, some combinations of initial solution and 
algorithm parameters may not lead to a moment-matching solution and the nonlinear 
programming algorithm may not converge. 

- For a given dimension we do not know how many solutions exist or how to guide 
the search toward a preferable solution. 

- The moment-matching solutions do not guarantee that the estimated cost function 
wUl have a finite minimum. When choosing among several solutions, a solution which 
gave a cost function with a finite minimum was always chosen. 

- The values of w, and W 2 can be modified to guide the search; generally w, ^ 1, 
W2^1 and w, ^ 3w2. The magnitude of w, and W 2 can be increased to get more or 
sufficient accuracy [Ref. 11: pp. 6]. 
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- If the feasible solutions are known, parameters can be bounded to guide the 
search. However, some bounds on the parameters may not lead to a moment-matching 
solution or to a solution whose cost function does not have a unique finite minim um. 

- User interaction is often necessary in obtaining the appropriate or preferable 
solutions. 

B. THE SEQUENTIAL ESTIMATION PROCEDURE 

Let Xi,X 2 , ,Xn be the sequence of lifetimes of the system and be the 

sequence of estimators of 0* where is the estimator after N replacements. Then after 
N observations, we have the right censored data (Zi,5i),(Z2,52),....,(Zn,6n)- 
where Z; = min (Xi,$*i_i) 

6i = 1 if Xj ^ $'i.i otherwise d, = 0 . 

Because the data are right censored the usual method of moments approach for estimation 
needs to be modifted. Rather than use sample moments calculated from an empirical 
distribution, we use moments from a nonparametric estimator of F based on the censored 
data. From the right censored data after N observations, the procedure to compute the 
estimators {$*;} is developed as follows. 

1. Use the right censored data sequence (Zi, 5 i),(Z 2 , 62 ), --jCZnjSn) to estimate 
moments with F = l-F by the product-limit estimator: 
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(3.4) 



Fit) 



n 



ilt it) 



( 



N-i 

N-i+1 



« 



(i) 



where Z^^^, Z(n> are the order statistics of Z,, Zj,..., Z^ and 6 ( 2 ),..., 6 (nj are 

ordered according to the ordering of Z(,), Z(n>. Then calculate probabilities Pj, 

associated with Z(o, by 

P( = F(Z(0 - f (Z(i)) . (3.5) 

Note, since F has discontinuities only at Z(j) where 6 (j) = 1 , that Pj = 0 when 6 (jj = 0. 
Then estimates of the k* moment are given by 



E[X^] = ( 3 - 6 ) 



The estimate of the 2"** standardized moment is 



A _ s/E[X^] -iE[X])^ 

^[X] 



and the estimate of the S'"* standardized moment is 



^ = i[X^] -2E[X]^X^] -i-2(E[X] gj 

iy/&[X^] - (^[X] )2)3 

2. Estimate parameters of PH distribution (a,T), by matching the second and third 
standardized moments using the nonlinear programming approach. 

When executing the nonlinear programming algorithm using the model described 
in the previous section, the det(T) 0 constraint is replaced by a constraint on the 
expected lifetime. The reason for doing this is that we could not formulate an algebraic 
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constraint equivalent to the constraint det(T) ^ 0 . The problem can be solved by a 
constraint on the expected lifetime, i.e. by taking 

oT'^e = S[X] (3.9) 

yv 

where E[X] is calculated in step 1 . This constraint will make the search for preferable 
solutions easier. The initial solution is chosen by the user. Some initial solutions may or 
may not lead to a moment matching solution. Typically, the initial vector a is 
(l/m,l/m,...,l/m)or (1,0,. ..,0) and an initial matrix T consists of tg = -1 and tjj = 1/m 
for i j, i = l,2,...,m. In this thesis the initial vector a is (1,0,..., 0) and an initial 
matrix T is taken to have tj j = -1 and tjj = 1/m [Ref. 11: pp. 9] . 

3. Using the parameters estimated in step 2, the estimated cost function is 



. ^ Cj - aexp(g^)e[C2-Cj 

[exp (t4>) - J] e 



(3.10) 



The new estimate of the optimal age of replacement is taken to be the <t> that 
minimizes (3.10) by enumerative search. 

4. Compare the with X^+i then repeat Step 1. 

A 

The initial solutions for matrix T and vector a in step 2 are taken to be the previously 
estimated values of T and a. In case the moment matching solution cannot be met, the 
initial solution in step 2 will be taken to be the original initial solution and the previous 
optimal age rqjlacement is used for step 3. 

The rq)lacement cost for i* system is Cj if X; < otherwise the replacement 
cost is Cl. The actual total rq>lacement cost for the first N systems that are observed is 
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(3.11) 



N 

= Yt [C'2x6_i + qx(i-6_i) ] 

1*1 

and the total operating time for the N systems is 

• (3.12) 

i-l i=l 

So, the actual average cost (AAC) after N replacement is computed by 

AAC^ = (3.13) 

In this thesis, the moment matching nonlinear programming approach for estimating 
parameters of PH distributions uses the GAMS program as shown in Appendix A (for 
a PH distribution with 3 transient states). The Fortran programs are written to estimate 
the optimal age of replacement, second and third standardized moments and for data 
preparation as shown in Appendbc B. 
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rV.COMPARISON WITH NONPARAMETRIC PROCEDURE 



A. NONPARAMETRIC PROCEDURE 

An alternate approach for estimating 0* is to estimate nonparametrically, as in Aras 
and Whitaker (1991). The procedure is much simpler computationally than the parametric 
procedure described in the previous section for PH distributions. In the nonparametric 

/V 

procedure after N replacements the product limit estimator F of F given in (3.4) is used 
to estimate the cost function C(</>) as 



f F{u)du 
J 0 



( 4 . 1 ) 



where F = 1 - F is the estimator of the cdf F, the estimator of 0* is then found by 
minimizing Cn(<^). In general one would expect parametric procedures to do much better 
than nonparametric procedures. However, with the numerical difficulties involved with 
estimating parameters for PH distribution and the fact that the family of PH distributions 

A. 

is so large it is not obvious which procedure is best. The criterion for comparison is the 
actual average cost per unit time, AAC^, after N rq>lacements. 



19 



B. COMPARISON OF NONPARAMETRIC AND PARAMETRIC PROCEDURE 



FOR THE PH DISTRIBUTION 

Simulation was used to compare sequential estimation based on PH distributions 
with nonparametric sequential estimation, Cj = 100, C 2 = 500. System lifetimes were 
generated by an increasing failure rate PH distribution with representation (T 4 i,a 4 i) 
where. 
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I - 6 5 

II - 7 



= (1 . 0 , 0 . 0 , 0 . 0 ) 



For this representation the long run expected cost per unit time is given in Figure 3. 

To give the parametric procedure a chance, the first 25 system lifetimes are 
uncensored. After the first 25 observations, the parameters of the PH distribution are 
estimated sequentially as in Chapter HI Section B. The actual average costs (AAC) are 
calculated and compared for small, medium and large sample sizes with planned cost Cj 
= 100 and unplanned cost C 2 = 500. The 20 initial data sets are simulated as shown in 
the Appendix A. 

The programs are separated into 5 parts. The first one is written in GAMS as 
shown in Appendix A, to estimate the parameters of a PH distribution by the nonlinear 
programming approach. The rest are written in Fortran as shown in Appendix B, to 
estimate the optimal age of replacement based on the estimated PH distribution from the 
GAMS program; to simulate the system lifetime; to prepare the right censored data; to 
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estimate the expected lifetime, second and third standardized moments; and finally to 
tabulate the results and to prepare the initial data for the next estimation. 



LONG RUN EXPECTED COST PER UNIT TIME OF PH DISTN. 
WITH REPRESENTATION T41, a41 AND 01 = 100, 02=500 




0 12 3 4 

REPUCEMENT AGE 



Figure 3 Long run expected cost per unit time curve of PH distribution with 
representation T 41 and 0 ( 4 , 

In this simulation, the number of transient states is fixed at 3 and the initial 
probability vector a = (1, 0, 0). By fixing a = (1, 0, 0) we reduce the number of 
parameters to be estimated and use a model in which aU systems start in the same state, 
when states are thought of as the level of system repair, this choice better reflects the 
idea that replacement systems are as good as new. Bounded parameters and user 
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interaction are used in guiding each estimation during the simulation to a preferable 
solution. 

Tables 1, 2 and 3 summarize the simulation results of the actual average cost 
resulting from sequential estimation of PH and the nonparametric procedure. Also given 
are the signed ranks of the differences in the actual average costs. These are used in the 
Wilcoxon Signed-Rank test to test the hypothesis: 

H„: E[AACp|J — E[AACnonp] 

H,: E[AACpjJ < E[AAC}|onp] 

by using T"^ the sum of the positive ranks as the test statistic . 

Tables 1, 2 and 3 give the statistic T"^ as 35, 10 and 2 respectively. When compared to 
the value in Table 9 [Ref. 13: pp. 780] at N=20 all p-values are less than 0.005 which 
implies that the sequential estimation of a PH distribution is better than sequential 
estimation of nonparametric for all sample sizes (small, medium and large). The box 
plots in Figures 5 and 6 of actual average cost versus the number of replacements at 
different sample sizes show that the actual average costs decrease as the number of 
replacements increase and that AACpH for the parametric procedure decrease more than 
the AACnonp for the nonparametric procedure. 

A second, PH distribution was chosen to compare the parametric PH procedure 
with the nonparametric procedure. The second PH distribution has an average long run 
expected cost function that is more shallow than the first PH distribution. It has 
representation T42 and «42 



22 



where 



^42 



-5 4 1 
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0^2 = ( 1 . 0 , 0 . 0 , 0 . 0 ) , 



and the long run expected cost per unit time is given in Figure 4. 



O 

to 

O 

o 

o 

Ui 

I— 

o 

UI 

OL O 

^ g 

Ui ^ 

2 

3 

q : 

O 

2 

S o 

tn 

in 



LONG RUN EXPECTED COST PER UNIT TIME OF PH DISTN. 
WITH REPRESENTATION T42, o42 AND Ct = 100. C2 = 500 
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Figure 4 Long run expected cost per unit time curve of PH distribution with 
representation T42 and 0^42 



Tables 4, 5 and 6 summarize the result and give the statistic of 93, 81 and 47 
respectively. The p-value taken from table 9 [Ref. 13: pp. 780] indicated that the 
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parametric procedure is not better than the nonparametric procedure for small and 
moderate sample sizes. In fact, for some cases (e.g. runs 16 and 20 in Table 4 and runs 
1, 5 and 16 in Table 5) the actual average cost resulting from using the nonparametric 
procedure can be quite a bit less than the actual average cost resulting from the 
parametric procedure. This is due to the fact that C(<^) is shallow at <l>*, thus even though 
C{<f>) may be a reasonable estimator of C(4>), the variance of its minimizing value 0* is 
higher than for the previous PH distribution. This means that <f>‘ is often underestimated. 
Because C(<^) increases so rapidly to the left of <f>*, underestimating <{>* can increase the 
actual average costs dramatically. Another difficulty is that the estimated C(</>) may be 
relatively flat around <t>*. This causes numerical difficulties with the nonlinear 
programming algorithm. For large sample sizes the parametric procedure does do better 
than the nonparametric procedure. Here C(<^) estimates C(<^) more closely and it is easier 
for the nonlinear programming algorithm to identify the correct solution. 

In Figures 5 and 6, the improvement in actual average cost with sample size is 
shown for both the parametric and nonparametric procedure for the PH distribution with 
representation T41 and a4,. As expected, both procedures improve (actual average costs 
decrease) with sample size. However, the actual average cost for the parametric 
procedure decreases faster than for the nonparametric procedure. Both procedures exhibit 
considerable variability. The solutions have got still mixed that we should do more study 
on another representation until the more precise solutions come up. 
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Table 1 COMPARISON OF ACTUAL AVERAGE COSTS FOR SMALL SAMPLE 



SIZES (N = 50) OF PH DISTRIBUTION WITH T 4 , AND « 4 , 



Run II AACpj^ 




AACrf 


DIFFERENCE 


RANK 


1 1 690.095 


754.555 


759.142 


-64.46 


16 


1 2 1 714.728 


775.712 


780.872 


-80.644 


14 


9 


700.342 


684.379 


688.578 


15.963 


4 * 


4 


666.014 


725.331 


738.987 


-59.317 


12 


5 


714.330 


757.528 


735.058 $ 


-43.198 


8 


8 


544.322 


621.505 


621.507 


-77.183 


14 


7 


721.505 


810.149 


810.149 


-88.644 


20 


8 


565.854 


617.458 


636.556 


-51.604 


11 


9 


676.511 


725.394 


702.676 $ 


-48.883 


14 


10 


698.572 


766.182 


809.629 


-67.61 


14 


11 


814.989 


767.586 


801.150 


47.403 


9 


12 


644.921 


667.553 


672.116 


-22.632 


S 


13 


761.516 


686.794 


679.892 ! 


74.722 


18 * 


16 


608.477 


601.215 


617.711 


7.262 


1 * 


15 


706.591 


767.782 


782.198 


-61.191 


14 


16 


629.294 


667.805 


682.719 


-38.511 


6 


17 


688.572 


755.475 


748.353 $ 


-66.903 


16 


16 


636.963 


624.132 


665.344 


12.831 


3 * 


19 


583.226 


625.161 


609.036 $ 


-41.935 


7 


20 


644.065 


654.691 


676.005 


-10.626 


2 



AACrf = Actual average cost replacing when failure occurs 
DIFFERENCE = AACpn - AACnqnp 
* positive rank 
T+ = 35 ; 
p-value < 0.005 

$ AACrf is less than only AACnonp 
! AACrf is less than AACpH and AACnqn? 
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Table 2 COMPARISON OF ACTUAL AVERAGE COSTS FOR MEDIUM 



SAMPLE SIZES (N = 100) OF PH DISTRIBUTION WITH T 41 AND 



1 Run 


AACpH 


AACnonp 


AACrf 


DIFFERENCE 


RANK 


1 


648.127 


723.170 


732.263 


-75.043 


18 


2 


690.971 


786.554 


789.767 


-95.583 


20 


3 


649.886 


644.710 


683.933 


5.176 


2 * 


4 


696.858 


747.640 


755.982 


-50.782 


11 


9 


690.195 


688.749 


682.840 ! 


1.446 


1 * 


6 


577.742 


607.414 


610.850 


-29.672 


5 


7 


686.414 


766.365 


762.245 $ 


-79.951 


IJ 


9 


550.106 


612.046 


638.286 


-61.94 


14 


9 


637.268 


696.249 


689.552 $ 


-58.981 


12 


IG 


684.253 


730.138 


786.920 


-45.885 


16 


11 


802.039 


790.698 


807.530 


11.341 




12 


615.087 


687.719 


695.174 


-72.632 


17 


13 


677.302 


666.119 


674.830 


11.183 


3 


14 


586.838 


628.687 


638.451 


-41.849 


8 


15 


726.213 


749.507 


737.792 $ 


-23.294 


5 


16 


626.605 


694.232 


667.462 $ 


-67.627 


16 


17 


671.908 


732.053 


758.912 


-60.145 


13 


18 


632.954 


667.156 


727.346 


-34.202 


7 


19 


635.989 


679.399 


667.141 $ 


-43.41 


9 


20 


590.881 


657.064 


701.735 


-66.183 


15 



AACrf = Actual average cost replacing when failure occurs 
DIFFERENCE = AACpH - AACnon? 

* positive rank 
T^ = 10 ; 
p-value < 0.005 

$ AACrf is less than only AACnonp 
! AACrf is less than AACrh and AACnonp 
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Table 3 COMPARISON OF ACTUAL AVERAGE COSTS FOR LARGE SAMPLE 



SIZES (N = 200) OF PH DISTRIBUTION WITH T 4 , AND « 4 , 



Run 


AACpH 


AACnonp 


AACrp 


DIFFERENCE 


RANK 


1 


622.560 


655.891 


655.865 $ 


-33.331 


7 


2 


658.640 


711.831 


708.966 $ 


-53.191 


11 


3 


623.823 


635.065 


671.611 


-11.242 


4 


4 


678.141 


749.866 


754.490 


-71.725 


14 


5 


674.548 


714.291 


707.926 $ 


-39.743 


9 


8 


585.435 


645.920 


647.927 


-60.485 


13 


7 


647.302 


729.159 


733.860 


-81.857 


20 


8 


581.711 


607.068 


649.386 


-25.357 


5 


5 


638.465 


678.078 


719.904 


-81.439 


19 


10 


655.709 


717.569 


766.020 


-61.86 


14 


11 


760.478 


779.362 


762.986 $ 


-18.884 


4 


12 


589.088 


649.545 


681.004 


-60.457 


IS 


13 


696.081 


679.816 


712.761 


16.265 


2 * 


18 


602.032 


643.703 


649.832 


-47.8 


10 


15 


685.977 


767.133 


755.331 $ 


-81.156 


18 


18 


646.752 


723.719 


663.489 $ 


-76.967 


IS 


17 


680.457 


713.146 


754.616 


-32.689 


6 


18 


674.300 


697.304 


740.394 


-23.004 


4 


19 


636.649 


674.845 


670.146$ 


-38.196 


8 


20 


588.468 


654.994 


670.194 


-66.526 


15 



AACrf = Actual average cost rq)lacing when failure occurs 
DIFFERENCE = AACpH - AACnonp 
* positive rank 
T-" = 2 ; 
p-value < 0.005 

$ AACrf is less than only AACnonp 
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Table 4 COMPARISON OF ACTUAL AVERAGE COSTS FOR SMALL SAMPLE 



SIZES (N = 50) OF PH DISTRIBUTION WITH T^^ AND «42 



Run 


AACpH 


AACj^onp 


AACrp 


DIFFERENCE 


RANK 


1 


618.463 


561.050 


545.871 ! 


57.413 


16 * 


2 


604.098 


684.564 


686.650 


-80.466 


19 


3 


674.353 


664.878 


637.614 ! 


9.475 


5 * 


4 


605.528 


559.966 


575.315 


45.562 


It * 


5 


580.054 


521.136 


535.683 


58.918 


17 * 


6 


547.534 


553.458 


556.857 


-5.924 


4 


7 


632.508 


625.711 


625.711 


6.797 


5 * 


9 


632.510 


647.089 


651.264 


-14.579 


S 


9 


490.099 


578.531 


566.695 $ 


-88.432 


10 


10 


629.150 


668.268 


670.786 


-39.118 


13 


11 


535.032 


536.771 


535.194 


-1.739 


S 


12 


605.358 


644.458 


645.689 


-39.100 


11 


13 


519.537 


552.440 


552.443 


-32.903 


10 


14 


532.469 


553.664 


569.868 


-21.195 


8 


15 


671.880 


672.507 


672.160$ 


-0.627 


S 


16 


645.125 


571.650 


576.024 


73.475 


18 * 


17 


495.932 


533.334 


539.149 


-37.402 


11 


18 


648.741 


677.909 


682.898 


-29.168 


8 


19 


577.580 


576.101 


600.152 


1.479 


2 * 


20 


560.629 


519.418 


532.841 


41.211 


14 » 



AACrf = Actual average cost replacing when failure occurs 
DIFFERENCE = AACpH - AACnqnp 



* positive rank 
T+ = 93 ; 



p-value > 0.05 

$ AACrp is less than only AACnqnp 
1 ^^.ACpp is less than and 
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Table 5 COMPARISON OF ACTUAL AVERAGE COSTS FOR MEDIUM 
SAMPLE SIZES (N = 100) OF PH DISTRIBUTION WITH T 42 AND 



Run 


AACpH 


AACnonp 


AACrf 


DIFFERENCE 


RANK 


2 


573.310 


548.751 


539.868 ! 


24.559 


13 * 


2 


588.501 


633.798 


623.661 $ 


-45.297 


16 


3 


618.891 


626.148 


608.227 ! 


-7.257 


2 


4 


621.602 


590.532 


601.378 


31.070 


14 * 


5 


557.189 


521.040 


535.590 


36.149 


15 * 


6 


549.717 


562.877 


569.251 


-13.160 


9 


7 


651.656 


629.631 


630.823 


22.025 


1C * 


8 


615.798 


624.193 


628.347 


-8.395 


4 


9 


475.757 


574.707 


567.688 $ 


-98.95 


20 


1(5 


602.762 


626.458 


628.035 


-23.696 


12 


11 


535.426 


537.073 


549.888 


-1.647 


1 


12 


580.480 


598.790 


590.995 $ 


-18.310 


7 


13 


525.770 


597.208 


597.210 


-71.438 


19 


14 


518.669 


540.978 


553.372 


-22.309 


11 


15 


634.978 


627.571 


633.019 


7.407 


3 ♦ 


16 


675.027 


616.239 


618.041 


58.788 


18 * 


17 


535.911 


585.289 


587.998 


-49.378 


17 


18 


617.194 


636.912 


645.975 


-19.718 


9 


19 


550.468 


564.214 


579.831 


-13.746 


6 


20 


539.256 


520.034 


531.189 


19.222 


8 * 



AACrf = Actual average cost replacing when failure occurs 
DIFFERENCE = AACpH - AACnonp 
* positive rank 
T+ = 81 ; 
p-value > 0.05 

$ AACrj, is less than only AACnon? 

! AACrf is less than AACpH and AACnqnp 
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Table 6 COMPARISON OF ACTUAL AVERAGE COSTS FOR LARGE SAMPLE 



SIZES (N = 200) OF PH DISTRIBUTION WITH T 42 AND 



Run 


AACpH 


AACnonp 


AACrf 


DIFFERENCE 


RANK 


1 


556.373 


550.131 


544.302 ! 


6.242 


3 * 


2 


574.206 


616.553 


621.391 


-42.347 


17 


3 


578.172 


591.186 


606.516 


-13.014 


7 


8 


594.573 


579.915 


588.547 


14.658 


8 * 


5 


580.261 


551.976 


560.463 


28.285 


14 ♦ 


6 


526.351 


558.604 


573.294 


-32.253 


15 


7 


632.238 


622.079 


624.275 


10.159 


8 ♦ 


8 


612.338 


621.511 


610.932 ! 


-9.173 


5 


5 


460.437 


567.069 


563.305 $ 


-106.632 


20 


10 


557.424 


599.181 


602.122 


-41.757 


1© 


11 


555.711 


539.303 


551.773 


16.408 


g ★ 


12 


560.534 


584.128 


589.078 


-23.594 


12 


13 


576.939 


621.856 


621.857 


-44.917 


18 


14 


568.925 


595.619 


599.305 


-26.694 


If 


15 


562.097 


583.057 


595.668 


-20.960 


11 


16 


584.778 


604.211 


606.979 


-19.433 


10 


17 


549.979 


607.061 


610.878 


-57.082 


19 


16 


631.719 


628.295 


639.159 


3.424 


1 * 


19 


578.172 


577.936 


607.344 


0.236 


1 * 


20 


538.683 


530.002 


570.516 


8.681 


4 



AACrf = Actual average cost replacing when failure occurs 
DIFFERENCE = AACpH - AACnonp 
* positive rank 
T^ = 47 ; 

0.01 < p-value < 0.025 
$ AACrf is less than only AACnonp 
! AACrf is less than AACpH and AACnonp 
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Figure 5 The box plot of actual average cost for different sample sizes of 
phase type sequential procedure for T 4 , and 0(41 
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PHASE TYPE SEQUENTIAL ESTIMATION PROCDURE 




100 150 200 

SAMPLE NUMBER 



ACTUAL AVERAGE COST 




Figure 6 The box plot of actual average cost for different sample sizes of 
nonparametric sequential procedure for T 4 | and 
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NONPARAMETRIC SEQUENTIAL ESTIMATION PROCEDURE 




V.CONCLUSIONS AND RECOMMENDATIONS 



In this thesis, the age replacement policy has been considered. A system is replaced 
at the time of failure or at a scheduled replacement age whichever comes first. The 
replacement system is assumed to be as good as new. The objective is to achieve the 
minimum long run expected maintenance cost per unit time. The system lifetimes are 
assumed to have PH distributions. Estimating the parameters of a PH distribution is 
difficult because it may have many parameters and its density function involves a matrix 
exponential. Moment matching with a nonlinear programming approach is used for 
estimating the parameters. This method does not guarantee that the estimated cost 
function will have a unique and finite minimum. In this thesis we restric our attention to 
the case with initial probability vector a = (1, 0, 0). Sequential estimation using the 
phase type procedure gives smaller average costs than the nonparametric procedure for 
both small and large sample sizes when the underlying long run average cost function has 
a very distinct minimum at <^*. When this cost function is shallow around <}>*, the 
parametric procedure does not do better than the nonparametric procedure for small 
samples. In fact, for cost functions that are shallow, it is better not to use a maintenance 
policy i.e. to take ^* = c» , until the sample sizes are large enough to give reasonable 
estimates of the underlying parameters. The reason for this is that for such cost 
functions, overestimating does not increase the long run average cost much, at the 
same time the decrease in the amount of censoring with over estimating 0* greatly 
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improves estimation of the underlying distribution function F. Conversely 
underestimating <}>* for these cost functions causes a drastic increase in long run average 
costs and increase censoring making estimation very difficult. 

Recommendation for the future research 
The foUowing forms a list of future research initiatives: 

-Examine the impact of the product-limit estimator in estimating standardized 
moments on the rest of procedure. 

-Do more experimentation on different representations of phase type distributions 
to get more precise conclusion. 

-Look for other policies using phase-sensitive estimate and modified replacement 
costs over time. 

-The phase-type sequential estimations we have been so far assume the underlying 
lifetimes are iid and after replacement the system is new. To be more realistic, we can 
modify the initial probability vector, increase transition rates between states, or both, 
over time. 

-Use phase-type sequential estimation when the underlying life distribution F comes 
from Gamma, Weibull etc.^ that have increasing failure rate and compare with the 
noiqiarametric procedure. Since we can use the denseness of phase type distribution 
property to approximate the set of distribution with support on [0, a) i.e., the set of 
lifetimes, we may get a better method when the distribution of lifetimes is not known. 
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APPENDIX A. 



$TITLE PH DISTRIBUTION ESHMATION PARAMETERS 
* 

* gams and dollar control options 

* (See Appendices B & C) 

OPTIONS 

WORK = lOOOCK), 

UMCOL = 0 , UMROW = 0 , SOLPRINT = OFF , DECIMALS = 2 
RESLIM = 100, ITERLIM =100000, OPTCR = 0.0 , SEED = 3141; 

* DEFINITIONS AND DATA 

* This program uses nonlinear programming to estimate the 3-transient states of phase 

* type distribution parameters. Matching the second and third standardized moments are 

* used. All adjustable data are prepared by the "FILE SCALAR2" and uses the 

* "$INCLUDE" statement bring to the program. 

SINCLUDE "FILE SCALAR2" 

* MODEI^ 

VARIABLE 

Z objective function value 
A first entry of the matrix 
B second entry of the matrix 
C third entry of the matrix 
D fourth entry of the matrix 
E fifth entry of the matrix 
F sixth entry of the matrix 
G seventh entry of the matrix 
H eighth entry of the matrix 
I ninth entry of the matrix 
ALl initial probability of being in state 1 
AL2 initial probability of being in state 2 
AL3 initial probability of being in state 3 ; 

EQUATION 

OBJ Define objective function 

EQl Set initial probability of being in state 1 (nonnegative) 

EQ2 Set initial probability of being in state 2(nonnegative) 
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EQ3 Set initial probability of being in state 3 (nonnegative) 

EQ4 Sum over all of initial prob. less than or equal 1.0 

EQ5 Set the first diagonal entry to be nonpositive 

EQ6 Set the second diagonal entry to be nonpositive 

EQ7 Set the third diagonal entry to be nonpositive 

EQ8 Set off diagonal entries in the first row to be positive 

EQ9 Set off diagonal entries in the first row to be positive 

EQIO Set sum off diagonal of the first row less than or equal first diagonal 

EQll Set off diagonal entries in the second row to be positive 

EQ12 Set off diagonal entries in the second row to be positive 

EQ13 Set sum off diagonal of the second row less than or equal second diagonal 

EQ14 Set off diagonal entries in the third row to be positive 

EQ15 Set off diagonal entries in the third row to be positive 

EQ16 Set sum off diagonal of the third row less than or equal third diagonal 

EQ17 Define the expected lifetime ; 

* MINIMIZE 
OBJ.. 

Z =E= 

*** 3*SQR(C - C(0)) *** 

3*P0WER((((2*AL1 *(POWER((-(F*H) +E*I),2) +(C*H-B*I)*(F*G-D*I) 

+ (-(C*E) +B*F)*(-(E*G) +D*H) 

+(-(F*H) +E*I)*(C*H-B*I) + (C*H-B*I)*(-(C*G) + A*I) +(-(C*E) +B*F)*(B*G-A*H) 

+ (-(F*H) +E*I)*(-(C*E) +B*F) +(C*H-B*I)*(C*D-A*F) + (-(C*E) +B*F)*(-(B*D) + 

A*E)) 

+ 2*AL2*((F*G-D*I)*(-(F*H)+E*I)+(-(C*G)+A*I)*(F*G-D*I)+(-(E*G)+D*H)* 
(C*D-A*F) 

+ (F*G-D*I)*(C*H-B*I) +POWER((-(C*G) + A*I),2) +(C*D-A*F)*(B*G-A*H) 

+(F*G-D*I)*(-(C*E)+B*F)+(-(C*G)+A*I)*(C*D-A*F)+(C*D-A*F)*(-(B*D)+A*E)) 

+ 2* AL3*((-(E*G) +D*H)*(-(F*H) +E*I) + (B*G-A*H)*(F*G-D*I) + (-(B*D) + A*E)* 
(-(E*G)+D*H) 

+ (-(E*G) +D*H)*(C*H-B*I) + (B*G-A*H)*(-(C*G) + A*I) +(-(B*D) + A*E)* 

(B*G-A*H) + (-(E*G) +D*H)*(-(C*E) +B*F) + (B*G-A*H)*(C*D-A*F) + 

POWER((-(B*D) + A*E),2)))/ 

POWER((-(C*E*G) + B*F*G + C*D*H - A*F*H - B*D*I + A*E*I),2) 



-POWER((-(ALl *(-F*H+ E*I+ C*H-B*I-(C*E) +B*F) 

+AL2*(F*G-D*I-C*G+ A*I+C*D-A*F)+AU*(-E*G+D*H+B*G-A’^-B*D+A*E))/ 
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(-(C-"E*G)+B-"F’^G+C-^D-^H-A*F-*‘H-B*D*I+A*E*I)),2))-"’^0.5/ 



(-(ALl *(-(F*H) +E*I+C*H-B*I-(C*E) +B*F) 

+AL2*(F*G-D*I-C*G+A*I+C*D-A*F) 

+ AL3*(-E*G +D“^H+B*G-A*H-B*D + A-^E))/ 

(-(C*E*G) + B*F*G + C*D*H - A*F*H - B*D*I + A*E*I)) - MM2), 2) 

*** W2*SQR(GAMMA - GAMMA(O)) *** 

+ P0WER(((((-6*AL1 *((POWER((-(F*H) +E*I) ,2) + (C’"H-B*I)*(F*G-D*I) + 

(-(C*E) +B-^F)*(-(E*G) +D»H))*(-(F»H) +E“^I) 

+ ((-(F’^ +E*I)*(C*H-B*I) + (C’^-B’^*(-(C*G) + A*^ + (-(C*E) +B*F)* (B*G-A*H))* 
(F^G-D-^I) 

+ ((-(F*H) +E*I) *(-(C*E) +B*F) + (C*H-B*I) -^(C^D-A^F) + (-(C*E) +B-^F)* 

(-(B*D) + A*E))*(-(E*G) + D*H) 

+ (POWER((-(F*H) +E*I),2) + (C*H-B*I)’^(F*G-D’^ + (-(C*E) +B*F)* (-(E*G) +D*H))* 
(C*H-B“^I) 

+ ((-(F*H) + E*I) *(C*H-B*I) + (C*H-B*I) ’X-(C*G) + A*I) + (-(C*E) +B*F)*(B*G-A-^H))* 
(-(C*G)+A*I) 

+ ((-(F-^H) +E*I)*(-(C*E) +B*F)+ (C*H-B*I)*(C*D-A*F) + (-(C*E) +B*F)* 

(-(B’^D) + A*E))*(B*G-A*H) 

+ (POWER((-(F*H) +E*I),2) +(C*H-B’"I)’"(F*G-D*I) +(-(C*E) +B*F)* 

(-(E*G) +D*H))*(-(C*E) +B*F) 

+ ((-(F*H) + E*I) *(C*H-B“"I) + (C*H-B*I)*(-(C*G) + A*I) + (-(C*E) +B*F)* 
(B*G-A’^H))’^(C*D-A*F) 

+ ((-(F*H) +E*I)*(-(C*E) +B*F) +(C*H-B*I)*(C*D-A*F) + (-(C*E) +B’"F)’^ 

(-(B*D) + A*E))*(-(B*D) + A*E)) 

-6*AL2’^(((F’^G-D*I)*(-(F*H) +E*I) + (-(C*G) + A*I) *(F*G-D*I) + (-(E*G) +D*H)* 
(C*D-A*F))*(-(F*H) +E’^I) 

+((F*G-D*I)*(C*H-B*I)+POWER((-(C*G)+A*I),2)+(C*D-A*F)*(B*G-A*H))* 

(F»G-D*I) 

+ ((F*G-D*I)*(-(C’"E) +B*F) + (-(C*G) + A*I)*(C*D-A*F) + (C*D-A*F)* 
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(-(B*D) + A*E))*(-(E*G) +D*H) 



+((F*G-D*I)*(-(F*H) +E*I) + (-(C*G) + A*I)*(F*G-D*I) +(-(E*G) +D*H)* 
(C*D-A*F))*(C*H-B*I) 

+((F*G-D*I)*(C*H-B*I)+P0WER((-(C*G)+A*I),2)+(C*D-A*F)*(B*G-A*H))* 

(-(C*G)+A*I) 

+ ((F*G-D*I)*(-(C*E) +B*F) + (-(C*G) + A*I)*(C*D-A*F) + (C*D-A*F)* 

(-(B*D) + A*E))*(B*G-A*H) 

+ ((F*G-D*I)*(-(F*H) +E*I) +(-(C*G) + A*I)*(F*G-D*I) + (-(E*G) +D*H)* 
(C*D-A*F))*(-(C*E) +B*F) 

+ ((F*G-D*I)*(C*H-B*I) +POWER((-(C*G) + A*I),2) + (C*D-A*F)*(B*G-A*H)) * 
(C*D-A*F) 

+ ((F*G-D*I)*(-(C*E) +B*F)+ (-(C*G) + A*I)*(C*D-A*F) + (C*D-A*F)* 

(-(B*D) + A*E))*(-(B*D) + A*E)) 

-6*AL3*(((-(E*G)+D*H)*(-(F*H)+E*I)+(B*G-A*H)*(F*G-D*I)+(-(B*D)+A*E)* 
(-(E*G) +D*H))*(-(F*H) +E*I) 

+ ((-(E*G) +D*H)*(C*H-B*I) + (B*G-A*H)*(-(C*G) + A*I) + (-(B*D) + A*E) * 
(B*G-A*H))*(F*G-D*I) 

+ ((-(E*G) +D*H)*(-(C*E) +B*F) + (B*G-A*H)*(C*D-A*F) + 

POWER((-(B*D) + A*E),2))*(-(E*G) +D*H) 

+ ((-(E*G) +D*H)*(-(F*H) +E*I) + (B*G-A*H)*(F*G-D*I) +(-(B*D) + A*E)* 
(-(E*G) +D*H))*(C*H-B*I) 

+ ((-(E*G) +D*H)*(C*H-B*I) +(B*G-A*H)*(-(C*G) + A*I) + (-(B*D) + A*E)» 
(B*G-A*H))*(-(C*G) + A*I) 

+ ((-(E*G) +D*H)*(-(C*E) +B*F) + (B*G-A*H)*(C*D-A*F) + 

POWER((-(B*D) + A*E),2))*(B*G-A*H) 

+ ((-(E*G) +D*H)*(-(F*H) +E*I) + (B*G-A*H)*(F*G-D*I)+ (-(B*D) + A*E)* 
(-(E*G)+D*H))*(-(C*E) +B*F) 

+ ((-(E*G) + D*H) *(C*H-B*I) + (B*G- A*H)*(-(C*G) + A*I) + (-(B*D) + A*E) * 
(B*G-A*H))*(C*D-A*F) 
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+((-(E*G) +D*H)*(-(C*E) +B*F) +(B*G-A*H)*(C*D-A*F) + 

POWER((-(B*D) + A*E) ,2))*(-(B*D) + A*E))) 

/POWER((-(C*E*G) + B*F*G + C*D*H - A*F*H - B*D*I + A’^E*I),3) 

+3*((AL1*(-(F*H)+E*I+C*H-B*I-(C*E)+B*F)+AL2*(F*G-D*I-C*G+A*I+C*D- 

A*F)+AL3*(-E*G+D*H+B*G-A*H-B*D+A*E))/(-C*E*G+B*F*G+C*D*H-A*F 

*H-B*D*I+A*E*I))* 



*** expected of lifetime square *** 

((2*AL1 *(POWER((-(F*H) +E*I),2) + (C*H-B*I)*(F*G-D*I) +(-(C*E) +B*F)* 
(-(E*G) +D*H) +(-(F*H) +E*I)*(C*H-B*I) +(C*H-B*I)*(-(C*G) + A*I) + 
(-(C*E)+B*F)*(B*G-A*H)+(-(F*H)+E*I)*(-(C*E)+B*F)+(C*H-B*I)*(C*D-A*F) 
+ (-(C*E) +B*F)*(-(B*D) + A*E)) 

+2*AL2*((F*G-D*I)*(-(F*H)+E*I)+(-(C*G)+A*I)*(F*G-D*I) 

+(-(E*G) +D*H)*(C*D-A*F)+(F*G-D*I)*(C*H-B*I)+POWER((-(C*G)+A*I),2) 
+(C*D-A*F)*(B*G-A*H)+(F*G-D*I)*(-(C*E)+B*F)+(-(C*G)+A*I)*(C*D-A*F) 
+(C*D-A*F)*(-(B*D)+A*E)) 

+2*AL3*((-(E*G) +D*H)*(-(F*H) +E*I) +(B*G-A*H)*(F*G-D*I) 

+(-(B*D) + A*E)*(-(E*G) +D*H) + (-(E*G) +D*H)*(C*H-B*I) +(B*G-A*H)* 

(-(C*G) + A*I) +(-(B*D) + A*E)*(B*G-A*H) + (-(E*G)+D*H)*(-(C*E) +B*F) 

+ (B*G-A*H)*(C*D-A*F) +POWER((-(B*D) + A*E),2)))/ 

POWER((-(C*E*G) + B*F*G + C*D*H - A*F*H - B*D*I + A*E*I),2)) 

*** 2 time cube of the expected lifetime *** 

-2*POWER(((ALl *(-(F*H) +E*I +C*H-B*I-(C*E) +B*F) 

+AL2*(F*G-D*I-C*G+A*I+C*D-A*F)+AL3*(-E*G+D*H+B*G-A*H-B*D+A*E))/ 
(-(C*E*G) + B*F*G + C*D*H - A*F*H - B*D*I + A*E*I)),3)) 

/((2*AL1 *(POWER((-(F*H) +E*I),2) + (C *H-B*I)*(F*G-D*I) + (-(C*E) +B*F)* 
(-(E*G) +D*H) +(-(F*H) +E*I)*(C*H-B*I) +(C*H-B*I)*(-(C*G) + A*I) + 
(-(C*E)+B*F)*(B*G-A*H)+(-(F*H)+E*I)*(-(C*E)+B*F)+(C*H-B*I)*(C*D-A*F) 
+(-(C*E) +B*F)*(-(B*D) + A*E)) 

+2*AL2*((F*G-D*I)*(-(F*H) +E*I) + (-(C*G) + A*I)*(F*G-D*I) + 

(-(E*G) +D*H)*(C*D-A*F) + (F*G-D*I) *(C*H-B*I) + 
POWER((-(C*G)+A*I),2)+(C*D-A*F)*(B*G-A*H)+(F*G-D*I)* 

(-(C*E) +B*F) + (-(C*G) + A*I)*(C*D-A*F) +(C*D-A*F)*(-(B*D) + A*E)) 
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+ 2*AL3*((-(E*G) +D*H)*(-(F*H) +E*I) + (B*G-A*H)*(F*G-D*I) + 

(-(B*D) + A*E)*(-(E*G) + D*H) + (-(E*G) + D*H)*(C*H-B*I) + (B*G-A’^H)* 
(-(C*G) + A*I) + (-(B*D) + A*E)*(B*G-A*H) + (-(E*G) +D*H)*(-(C*E) +B*F) 
+(B*G-A*H)*(C*D-A*F) +POWER((-(B*D) + A*E),2)))/ 

(POWER((-(C*E*G) + B*F*G + C*D*H - A*F*H - B*D*I + A*E*I),2)) 

-P0WER(((AL1*(-(F*H) +E*I + C*H-B*I-(C*E) +B*F) 
+AL2*(F*G-D*I-C*G+A*I+C*D-A*F) 

+ AL3*(-E*G +D*H +B*G-A*H-B*D + A*E))/ 

(-(C*E*G) + B*F*G + C*D*H - A*F*H - B*D*I + A*E*I)),2))**1.5) - MM3), 2); 
♦ 

* SUBJECT TO 

♦ 

EQl.. 

ALl =G= 0 ; 

EQ2.. 

AL2 =G= 0 ; 

EQ3.. 

AL3 =G= 0 ; 

EQ4.. 

ALl + AL2 + AL3 =L= 1 ; 

EQ5.. 

A =L= 0 ; 

EQ6.. 

E =L= 0 ; 

EQ7.. 

I =L= 0 ; 

EQ8.. 

B =G= 0 ; 

EQ9.. 

C =G= 0 ; 

EQIO.. 

B + C =L= -A ; 

EQll.. 

D =G= 0 ; 

EQ12.. 

F =G= 0 ; 

EQ13.. 

D + F =L= -E ; 

EQ14.. 

G =G=0; 
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EQ15.. 

H =G=0; 

EQ16.. 

G + H =L= -I ; 

EQ17.. 

-(ALl *(-(F*H) +E*I+ C-^H-B*I-(C*E) -f B*F) 
+AL2*(F*G-D*I-C*G+A*H-C*D-A*F) 
+AL3*(-E’^G+D*H+B*G-A’^H-B*D+A*E))/ 
(-(C*E*G)+B*F*G+C*D*H-A*F*H-B*D*I+A*E*I) =E= EX ; 

♦ 

* REPORT- 

MODEL MATCHM /ALL/ ; 

* Parameters bounding 

A.LO = -7.0; B.LO = 3.0; C.LO = 0.0; D.LO = 0.0; E.LO =-7.0; 

A.UP = -5.0; B.UP = 5.0; C.UP = 2.0; D.UP = 1.0; E.UP =-5.0; 

F.LO = 4.0; G.LO = 0.0; H.LO = 0.0; I.LO = -8.0; 

F.UP = 5.0; G.UP = 1.0; H.UP = 2.0; I.UP = -7.0; 

ALl.LO = 1.0; AL2.LO = 0.0; AL3.LO = 0.0; 

ALl.UP = 1.0; AL2.UP = 0.0; AL3.UP = 0.0; 

* The initial solution 

A.L = AA ; B.L = BB ; C.L = CC; D.L = DD; E.L = EE ; 

F.L = FF ; G.L = GG ; H.L = HH ; I.L = H ; 

ALl.L = ALPl ; AL2.L = ALP2 ; AL3.L = ALP3 ; 

SOLVE MATCHM USING NLP MINIMIZING Z ; 

DISPLAY Z.L,A.L , B.L, C.L, D.L, E.L, F.L, G.L, H.L, I.L, 
ALl.L, AL2.L, AL3.L; 

* Write output to the CMS file 

FILE OUTl /FILE OBJ A / ; 

PUT OUTl ; 

PUT Z.L ; 

FILE OUT2 /FILE ALPHA A / ; 

PUTOUT2 ; 

PUT A.L ; 

FILE OUT3 /FILE BRAVO A / ; 

PUTOUT3 ; 
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PUT B.L ; 

FILE 0UT4 /FILE CHAREE A / ; 
PUT 0UT4 ; 

PUT C.L ; 

FILE OUTS /FILE DELTA A / ; 
PUT OUTS ; 

PUT D.L ; 

FILE OUT6 /FILE ECHO A / ; 
PUT OUT6 ; 

PUT E.L ; 

FILE OUT7 /FILE FOXTROT A / ; 
PUT OUT7 ; 

PUT F.L ; 

FILE OUT8 /FILE GOLF A / ; 

PUT OUTS ; 

PUT G.L ; 

FILE OUT9 /FILE HOTEL A / ; 
PUTOUT9 ; 

PUT H.L ; 

FILE OUTIO /FILE INDIA A / ; 
PUT OUTIO ; 

PUT I.L ; 

FILE OUTll /FILE ALl A / ; 

PUT OUTll ; 

PUT ALl.L ; 

FILE OUT12 /FILE AL2 A / ; 

PUT OUT12 ; 

PUT AL2.L ; 

FILE OUT13 /FILE AL3 A / ; 

PUT OUT13 ; 

PUT AL3.L ; 
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THE INITIAL COMPLETE SYSTEM LIFETIME SET 1-5 



N 


SET 1 


SET 2 


SET 3 


SET 4 


SET 5 


1 


0.1226 


0.0740 


0.6855 


0.0528 


0.1143 


2 


0.1455 


0.1413 


0.2157 


0.1139 


0.1999 


3 


0.2662 


0.2261 


0.2209 


0.2412 


0.2002 


a 


0.2925 


0.2629 


0.2839 


0.2892 


0.2422 


5 


0.3205 


0.2831 


0.2852 


0.3118 


0.3164 


5 


0.4214 


0.3340 


0.3116 


0.3485 


0.3398 


7 


0.4394 


0.3911 


0.3403 


0.3973 


0.3534 


8 


0.4685 


0.4517 


0.4517 


0.4442 


0.3893 


8 


0.5804 


0.4761 


0.4691 


0.4686 


0.4410 


10 


0.5871 


0.4909 


0.5254 


0.4777 


0.4766 


11 


0.6228 


0.5290 


0.5536 


0.4885 


0.4964 


12 


0.6534 


0.5818 


0.5993 


0.4970 


0.5074 


15 


0.7220 


0.5962 


0.6575 


0.5114 


0.5240 


10 


0.7590 


0.6233 


0.6623 


0.6073 


0.5825 


If 


0.7618 


0.6419 


0.6855 


0.6433 


0.6033 


16 


0.8206 


0.7407 


0.7156 


0.7576 


0.6255 


17 


0.8245 


0.7682 


0.7690 


0.8057 


0.7391 


le 


0.8317 


0.8312 


0.7866 


0.8740 


0.7723 


19 


1.0056 


0.8693 


0.2852 


1.IL004 


0.9128 


20 


1.0153 


0.9727 


0.9473 


1.0424 


1.0832 


21 


1.0853 


1.1084 


1.0478 


1.0634 


1.1936 


22 


1.1659 


1.1272 


1.1155 


1.0659 


1.2003 


23 


1.1852 


1.2037 


1.2936 


1.2238 


1.8976 


24 


1.2040 


1.3319 


1.3906 


1.0424 


2.0769 


25 


1.3069 


2.5362 


1.6715 


3.0033 


2.2631 1 


E[X] 


0.705 


0.704 


0.679 


0.735 


0.742 1 


2“* MM 


0.47917 


0.70936 


0.57558 


0.80566 


0.77290 1 


3"* MM 


0.03190 


1.88276 


0.73612 


2.23985 


1.40912 I 
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THE INITIAL COMPLETE SYSTEM LIFETIME SET 6-10 



N 


SET 6 


SET 7 


SET 8 


SET 9 


SET 10 


1 


0.1337 


0.0402 


0.1759 


0.0830 


0.0761 


2 


0.1538 


0.0634 


0.2031 


0.1101 


0.0913 


3 


0.2452 


0.1939 


0.2106 


0.2083 


0.2235 


4 


0.2876 


0.1954 


0.2376 


0.2339 


0.3309 


9 


0.3591 


0.1958 


0.3043 


0.2408 


0.3470 


8 


0.4589 


0.2353 


0.3719 


0.2726 


0.3662 


7 


0.5161 


0.2775 


0.3850 


0.2777 


0.4091 


8 


0.6030 


0.3185 


0.3892 


0.2894 


0.4523 


9 


0.6181 


0.3579 


0.4032 


0.3137 


0.5581 


13 


0.6761 


0.3922 


0.4199 


0.3227 


0.5800 


18 


0.6851 


0.3985 


0.6029 


0.3390 


0.5805 


12 


0.7738 


0.5206 


0.6406 


0.387t 


0.5993 


13 


0.8541 


0.5233 


0.6506 


0.3878 


0.6354 


18 


0.8654 


0.6261 


0.7156 


0.4055 


0.6571 


18 


0.4589 


0.6336 


0.7584 


0.4076 


0.7529 


18 


1.0193 


0.6776 


0.7677 


0.4333 


0.7874 


17 


1.0352 


0.8705 


0.8754 


0.6178 


0.7991 


18 


1.0677 


1.1124 


0.6836 


0.6245 


0.8572 


18 


1.1196 


1.1228 


0.9575 


0.6455 


0.8920 


20 


1.1995 


1.1742 


0.9805 


0.9824 


0.9024 


21 


1.2301 


1.1872 


1.3367 


1.0054 


0.9699 


22 


1.3043 


1.1887 


1.4040 


1.1143 


0.9734 


23 


2.0273 


1.4849 


1.5844 


1.2279 


1.0452 


20 


3.2891 


1.7218 


1.6250 


1.2420 


1.1461 


25 


3.5360 


2.5301 


1.8392 


2.1319 


1.1594 


E[X] 


1.000 


0.722 


0.749 


0.572 


0.648 


2““ MM 


0.82457 


0.81263 


0.63261 


0.80948 


0.46844 


3'’* MM 


1.89734 


1.25687 


0.79102 


1.69217 


0.14707 
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THE INITIAL COMPLETE SYSTEM LIFETIME SET 11-15 



I N 


SET 11 


SET 12 


SET 13 


SET 14 


SET 15 


1 


0.0268 


0.9433 


0.0470 


0.1730 


0.1101 


2 


0.0745 


0.1285 


0.0922 


0.2556 


0.1386 


9 


0.1235 


0.2052 


0.1295 


0.3124 


0.1709 


1 


0.1800 


0.2435 


0.2064 


0.3125 


0.1998 


1 


0.2526 


0.2924 


0.2924 


0.3447 


0.2075 


8 


0.2704 


0.3335 


0.2957 


0.3894 


0.2588 


t 


0.2972 


0.3645 


0.3110 


0.4044 


0.2669 


8 


0.3697 


0.4936 


0.3209 


0.5559 


0.2737 


9 


0.3811 


0.4984 


0.4173 


0.5796 


0.2501 


10 


0.3990 


0.5064 


0.4577 


0.5803 


0.3506 


11 


0.4070 


0.5778 


0.5168 


0.5928 


0.3604 


12 


0.4855 


0.6699 


0.5179 


0.6056 


0.4117 


13 


0.6402 


0.7442 


0.5258 


1.6525 


0.5809 


10 


0.7087 


0.7502 


0.6056 


0.7299 


0.5981 


15 


0.7408 


0.8656 


0.5547 


0.7355 


0.6558 


18 


0.8496 


0.8831 


0.5559 


0.8122 


0.7306 


It 


0.9290 


1.0607 


0.5633 


0.9379 


0.7388 


12 


0.9873 


1.1003 


0.5692 


0.9433 


0.7400 


15 


0.9983 


1.1834 


0.6501 


1.0197 


0.7786 


20 


1.0318 


1.2235 


0.7978 


1.2675 


0.7936 


21 


1.0937 


1.3370 


0.8750 


1.3842 


1.1532 


22 


1.1137 


1.5754 


1.0843 


1.6525 


1.2444 


23 


1.2706 


1.8085 


1.1683 


1.8276 


1.2781 


24 


1.4875 


1.9023 


1.1881 


1.8545 


1.3428 


25 


1.4958 


1.9894 


1.3042 


2.8764 


1.5937 


E[X] 


0.665 


0.831 


0.560 


0.872 


0.611 


2“'* MM 


0.64986 


0.66778 


0.60540 


0.71398 


0.68553 


3"* MM ] 


0.32858 


0.60712 


0.69085 


1.51250 


0.77568 1 
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THE INITIAL COMPLETE SYSTEM LIFETIME SET 16-20 



N 


SET 16 


SET 17 


SET 18 


SET 19 


SET 20 


1 


0.1677 


0.0565 


0.0398 


0.1798 


0.0202 


2 


0.1819 


0.1811 


0.1402 


0.2037 


0.1920 


3 


0.2139 


0.1941 


0.1506 


0.2197 


0.2205 


3 


0.3436 


0.2810 


0.179t 


0.3159 


0.2209 


5 


0.3641 


0.3315 


0.1922 


0.3464 


0.3406 


6 


0.4393 


0.4147 


0.3787 


0.3629 


0.3645 


3 


0.4800 


0.4249 


0.6660 


0.3744 


0.3941 


9 


0.5328 


0.4629 


0.4944 


0.3793 


0.4130 


9 


0.5643 


0.9939 


0.5952 


0.3870 


0.5525 


IE 


0.7220 


0.5382 


0.6198 


0.3875 


0.5729 


11 


0.7393 


0.5776 


0.6269 


0.1011 


0.5787 


12 


0.7666 


0.6003 


0.6276 


0.4906 


0.5835 


13 


0.7766 


0.6055 


0.6608 


0.7436 


0.6138 


14 


0.7964 


0.6230 


0.7183 


0.7690 


0.6389 


15 


0.7981 


0.6477 


0.7970 


0.7824 


0.6609 


14 


0.^229 


0.7631 


0.7997 


0.8226 


0.6624 


17 


0.8390 


0.79?0 


0.8065 


0.8492 


0.7091 


11 


0.9422 


0.7988 


0.8065 


0.9542 


0.7178 


19 


1.0115 


0.8518 


0.9682 


1.2333 


0.7447 


20 


1.0381 


0.9091 


0.9939 


1.3601 


1.2027 


21 


1.1195 


1.0447 


1.0028 


1.4489 


1.3247 


22 


1.2716 


1.1973 


1.0875 


1.5989 


1.4412 


23 


1.3426 


1.2331 


1.2722 


2.2702 


1.6155 


24 


1.6082 


1.2339 


1.4913 


2.2971 


1.7162 


25 


1.7724 


1.3386 


1.8663 


3.1519 


2.1912 


E[X] 


0.786 


0.664 


0.711 


0.893 


0.748 


2““ MM 


0.52328 


0.51958 


0.60066 


0.83695 


0.70314 


3"* MM 


0.55452 


0.32242 


0.66254 


1.46868 


1.13073 1 
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APPENDIX B. 

PROGRAM COST 
C 

C THIS PROGRAM PROVIDES THE OPTIMAL AGE REPLACEMENT OF THE 
C PHASE TYPE DISTRIBUTION. IT WORKS ONLY SPECIHC CASE OF 3 
C PHASES OF TRANSIENT STATES. THE PARAMETERS ARE AS FOLLOWS: 
C T(I,J) = ENTRIES OF THE TRANSITION MATRIX OF PH DISTN. 

C AL(I) = INITIAL PROBABHJTY OF BEING IN STATE I 
C Cl = PLANNED MAINTENANCE COST 
C C2 = UNPLANNED MAINTENANCE COST 

C L = MAXIMUM LIFETIME THAT WANT TO CALCULATE 

C DET = DETERMINANT OF TRANSITION MATRIX 
C SMALL = THE OPTIMAL AGE REPLACEMENT 
C CS = AGE REPLACEMENT COSTS 

>K 

* INITIAL INPUT 

♦ 

INTEGER I,J,K,N 

REAL Cl, C2, TT(3,3), T(3,3), DET, AL(3), L, OBJ 
REAL TS(3,3), SMALL, PI, A, IDEN(3,3), S, NUMER 
REAL DENOM, P(3,3), Q(3,3), CS(400), EXPTS(3,3) 

DATA C1,C2,L,TS/100.0,500.0,4.0,9*0.0 / 

N = 0 
S = 0.01 

OPEN(UNIT =11, FILE = ’ALPHA’) 

OPEN(UNIT =12,FILE = ’BRAVO’) 

OPEN(UNIT =13,FILE = ’CHAREE’) 

OPEN(UNIT =14,FILE = ’DELTA’) 

OPEN(UNIT =15,FILE = ’ECHO’) 

OPEN(UNIT =16,FILE = ’FOXTROT’) 

OPEN(UNIT =17,FILE = ’GOLF’) 

OPEN(UNIT =18,FILE = ’HOTEL’) 

OPEN(UNIT =19,FILE = ’INDIA’) 

OPEN(UNIT =20,FILE = ’ALl’) 

OPEN(UNIT =21, FILE = ’AL2’) 

OPEN(UNIT =22,FILE = ’AL3’) 

OPEN(UNIT =23, FILE = ’OBJ’) 

READ(11,*)T(1,1) 

READ(12,*)T(1,2) 

READ(13,’^)T(1,3) 
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READ(14,*)T(2,1) 

READ(15,*)T(2,2) 

READ(16,*)T(2,3) 

READ(17,*)T(3,1) 

READ(18,*)T(3,2) 

READ(19,*)T(3,3) 

READ(20,*)AL(1) 

READ(21,*)AL(2) 

READ(22,*)AL(3) 

READ(23,*)OBJ 

IF(OBJ .GT. 0.05) GOTO 555 

♦ 

* CALCULATION COST FOR EACH REPLACEMENT TIME 

♦ 

DET = -(T(1,3)*T(2,2)*T(3,1)) 

/ + T(1,2)*T(2,3)*T(3,1) 

/ + T(1,3)*T(2,1)*T(3,2) 

/ -T(1,1)*T(2,3)*T(3,2) 

/ -T(1,2)*T(2,1)*T(3,3) 

/ + T(1,1)*T(2,2)*T(3,3) 

PRINT*,’DET = ’,DET 

TI(1,1) = (-(T(2,3)*T(3,2)) + T(2,2)*T(3,3))/DET 
TI(1,2) = (T(1,3)*T(3,2) - T(1,2)*T(3,3))/DET 
TI(1,3) = (-(T(1,3)*T(2,2)) + T(1,2)*T(2,3))/DET 
TI(2,1) = (T(2,3)*T(3,1)-T(2,1)*T(3,3))/DET 
TI(2,2) = (-(T(1,3)*T(3,1)) + T(1,1)*T(3,3))/DET 
Tl(2,3) = (T(1,3)*T(2,1)-T(1,1)*T(2,3))/DET 
TI(3,1) = (-(T(2,2)*T(3,1)) + T(2,1)*T(3,2))/DET 
11(3,2) = (T(1,2)*T(3,1)-T(1,1)*T(3,2))/DET 
11(3,3) = (-(T(1,2)*T(2,1)) + T(1,1)*T(2,2))/DET 
DO 10 I = 1,3 
DO 5 J = 1,3 

IF(I.EQ.J)THEN 
IDENa,J) = 1.0 
ELSE 

IDEN(I,J) = 0.0 
ENDIF 

5 CONTINUE 
10 CONTINUE 

♦ 

100 CONTINUE 
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N = N+1 
DO 20 I = 1,3 
DO 15 J = 1,3 

TS(I,J) = S*T(I,J) 

15 CONTINUE 

20 CONTINUE 

♦ 

* FIND THE EXPONENTIAL OF MATRIX TS 

* 

CALL EXPON(TS,EXPTS,PI) 

A = 0.0 
DO 21 I = 1,3 
DO 22 J = 1,3 

A = A + AL(I)*EXPTS(I,J) 

22 CONTINUE 

21 CONTINUE 

NUMER = C2 - A*(C2-C1) 

DO 301 = 1,3 
DO 25 J = 1,3 

P(I,J) = EXPTS(I,J) - IDEN(I,J) 

25 CONTINUE 
30 CONTINUE 
DO 45 I = 1,3 
DO40 J = 1,3 
Q(I,J) = 0.0 
DO 35 K = 1,3 

Q(I,J) = Qa,J) +na,K)*P(K,J) 

35 CONTINUE 
40 CONTINUE 
45 CONTINUE 
DENOM = 0.0 
DO 55 I = 1,3 
DO50 J = 1,3 

DENOM = DENOM + AL(I)*Q(I,J) 
50 CONTINUE 
55 CONTINUE 

CS(N) = NUMER/DENOM 
S = S+0.01 
IF(S.LE.L)GOTO 100 

* 

* FIND THE TIME THAT HAS MINIMUM COST 
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u u u u u 



K = 1 

200 IF(K.EQ.N)THEN 
SMALL = K-^0.01 
ELSEIF(CS(K+ 1).GT.CS(K))THEN 
SMALL = K*0.01 
ELSEIF(CS(K+ 1).LE.CS(K))TEIEN 
SMALL = (K+l)*0.01 
K = K + 1 
GOTO 200 
ENDIF 

♦ 

* OUTPUT OPTIMAL REPLACEMENT AGE 

* 

OPEN(UNIT = 3, FILE = ’OPTAGE’) 

WRITE(3,333)SMALL , PI 
333 F0RMAT(1X,F6.3,4X,E11.4) 

OPEN(UNIT = 4, FILE = ’OPTCOST’) 

WRITE(4,*)CS(K) 

555 STOP 
END 

SUBROUTINE EXPON(A,EXPA,PI) 

TfflS SUBROUTINE USES TO FIND THE EXPONENTIAL OF A 3X3 MATRIX. 
IT WORK WITH COUPLE OF IMSL SUBROUTINES 
(’EVCRG’ , ’UNCG’ , ’EPIRG’). 

INTEGER LDA,LDEVEC,N,NOUT,I,J,K 

PARAMETER(N = 3,LDA =N,LDEVEC =N,LDUINV =N,LD =N,LDU =N) 
REAL A(LDA,N),PI,EXPA(3,3) 

COMPLEX EVAL(N), EVEC(LDEVEC,N), ELD(3,3), UINV(3,3), 
COMPLEX C(3,3), D(3,3) 

OPEN(UNIT= 1 ,FILE= ’EXPONENT’) 

CALCULATE THE EIGENVALUES AND EIGENVECTERS 

« 

CALL EVCRG(N,A,LDA,EVAL,EVEC,LDEVEC) 

DO 15 I = 1,3 
DO 10 J = 1,3 
ELD(I,J) = 0.0 
10 CONTINUE 
15 CONTINUE 
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ELD(1,1) = CEXP(EVAL(1)) 

ELD(2,2) = CEXP(EVAL(2)) 

ELD(3,3) = CEXP(EVAL(3)) 

* 

* CALCULATE THE INVERSE EIGENVECTERS’ MATRIX 
■« 

CALL UNCG(N,EVEC,LDU,UINV,LDUINV) 

DO 30 I = 1,3 
DO 25 J = 1,3 

Ca,J) = 0.0 

DO20K = 1,3 

Ca,J) = Ca,J) + ELD(I,K)*UINV(K,J) 

20 CONTINUE 

25 CONTINUE 
30 CONTINUE 
DO 45 I = 1,3 
DO40 J = 1,3 
Da,J) = 0.0 

DO 35 K = 1,3 

D(I,J) = D(I,J) + EVEC(I,K)*C(K,J) 

35 CONTINUE 

40 CONTINUE 
45 CONTINUE 
DO 55 I = 1,3 
DO 50 J = 1,3 

EXPA(I,J) = REAL(D(I,J)) 

50 CONTINUE 
55 CONTINUE 
C 

C CALCULATE THE PERFORMANCE INDEX OF EIGENVALUES AND 
C EIGENVECTERS IF LESS THAN 1 ’EXCELLENT’, IF IT IS BETWEEN 1 
C AND 100 ’GOOD’, IF ITIS GREATER THAN 100 ’THE SOLUTION IS 
C SUSPECTED’ 

C 

PI = EPIRG(N,N,A,LDA,EVAL,EVEC,LDEVEC) 

RETURN 

END 
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PROGRAM SIMLIF 



C 

C TfflS PROGRAM PROVIDES A RANDOM NUMBER OF A PHASE TYPE 
C DISTRIBUTION THAT HAVE SPECIFIC REPRESENTATION WITH AN 
C INITIAL PROBABUJIIES VECTOR(ALPHA) AND A MATRIX OF 
C TRANSITION RATES BETWEEN TRANSIENT STATES. THIS PROGRAMS 
C WORKS WITH SUBROUTINE ’RANNUM’ AND MORE SPECIFIC CASE 
C ONLY 4 STATES DISTRIBUTION. THE VARIABLES ARE AS FOLLOWS: 

C LD(I,J) = TRANSITION RATE FROM STATE I TO STATE J 

C P(I,J) = TRANSITION PROBABILITY FROM STATE I TO STATE J 
C AL(I) = INITIAL PROBABILITY OF BEING IN STATE I 

C EUF = THE SUPPOSED UPPER BOUND EXPECTED LIFETIME 

* 

* INmAUZATTON 

* 

INTEGER S,SEED,N,N1 
REAL U,EXP,UF,ELIF 

REALLD12,LD13,LD14,LD21,LD23,LD24,LD31,LD32,LD34 
REAL P12,P13,P21,P23,P31,P32,AL1,AL2,AL3,D 
DATA LD12,LD13,LD14,LD21 ,LD23,LD24,LD31 ,LD32,LD34 
+ 75.0,5.0, 5.0, 6.0, 6.0, 6.0, 7.0, 7.0, 7.0/ 

DATA P12,P13,P21 ,P23,P31 ,P32,AL1 ,AL2,AL3,ELIF 
+ 70.8,0.2,0.167,0.833,0.143,0.286,1.0,0.0,0.0,10.0/ 

CALL EXCMS(’FILEDEF 10 DISK FILE LFTEST (DISP MOD’) 
OPEN(UNIT = 1,FILE = ’LIFETIME’) 

OPEN(UNIT = 2,FILE =’SEED’) 

READ(2,*)SEED 

CLOSE(2) 

LIF = 0.0 



* GENERATE UNIFORM 0 AND 1 FOR INITIAL PROBABILITY 

♦ 

CALL RANNUM(1,SEED,0.0,1.0,0,U) 

IF (U.LE.AL1) THEN 
S = 1 

ELSEIF(U.LE.AL1 + AL2)THEN 
S = 2 
ELSE 
S = 3 
ENDIF 

♦ 

* SIMULATE LIFETIME WHEN IS IN STATE 1 
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5 CONTINUE 

IF ((S.EQ.l) .AND. (LIF.LT.ELIF)) THEN 
CALL RANNUM(1, SEED, 0.0, 1.0,0,U) 
IF(U.GT.P12+P13)THEN 
S = 4 

CALL RANNUM(3,SEED,LD14,0,0,EXP) 
LIF = LIF + EXP 
ELSEIF(U.GT.P12)THEN 
S = 3 

CALL RANNUM(3,SEED,LD13,0,0,EXP) 
LIF = LIF + EXP 
ELSE 
S = 2 

CALL RANNUM(3,SEED,LD12,0,0,EXP) 
LIF = LIF + EXP 
ENDIF 

ENDIF 

* 

* SIMULATE LIFETIME WHEN IS IN STATE 2 

♦ 

IF((S.EQ.2) .AND. (LIF.LT.ELIF))THEN 
CALL RANNUM(1,SEED,0.0,1.0,0,U) 
IF(U.GT.P21 +P23)THEN 
S = 4 

CALL RANNUM(3,SEED,LD24,0,0,EXP) 
UF = LIF + EXP 
ELSEff(U.GT.P21)THEN 
S = 3 

CALL RANNUM(3,SEED,LD23,0,0,EXP) 
UF = UF + EXP 
ELSE 
S = 1 

CALL RANNUM(3,SEED,LD21,0,0,EXP) 
UF = UF + EXP 
ENDIF 

ENDIF 

K< 

* SIMULATION LIFETIME WHEN IS IN STATE 3 

♦ 

IF((S.EQ.3) .AND. (UF.LT.EUF))THEN 
CALL RANNUM(1,SEED,0.0,1.0,0,U) 
IF(U.GT.P31 +P32)THEN 
S = 4 
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CALL RANNUM(3, SEED, LD34, 0,0, EXP) 

LEF = LEF + EXP 
ELSEIF(U.GT.P3 1)TEIEN 
S = 2 

CALL RANNUM(3,SEED,LD32,0,0,EXP) 

LEF = LEF + EXP 
ELSE 
S = 1 

CALL RANNUM(3,SEED,LD31,0,0,EXP) 

LEF = LEF + EXP 
ENDIF 
ENDIF 

IF((S.LT.4) .AND. (LIF.LT.ELIF))GOTO 5 

♦ 

* OUTPUT TO ’FILE LIFETIME A’ 

♦ 

WRITE(10,*)LIF 

WRITE(1,*)LIF 

OPEN(UNIT = 2,FILE =’SEED’) 

WRITE(2,*)SEED 

STOP 

END 

SUBROUTINE RANNUM(DISTN,SEED,RPARM1,RPARM2,IPARM,X) 

C 

C THIS SUBROUTINE PROVIDES A RANDOM NUMBER OF 7 DISTRIBUTION. 
C IT INTERFACES WITH THE LLRANDOMH ROUTINES PROVIDED IN THE 
C NONIMSL LIBRARY. THE PARAMETER AND CALLING PROCEDURE ARE 
C AS FOLLOWS: 

C DIST = DISTRIBUTION TYPE YOU WANT TO SELECT AN INTEGER 
C BETWEEN 1 AND 7. 

C SEED = THE RANDOM NUMBER SEED YOU WISH TO USE. 

C RPARMl , RPARM2, AND IPARM ARE REAL AND INTEGER PARAMETERS. 
C PASSED TO THE ROUTINE WITH MEANINGS WHICH VARY WITH THE 
C TYPE OF DISTRIBUTION YOU DESIRE. 

C X = THE RETURNED RANDOM NUMBER, IT IS ALWAYS REAL. 

C DISTRIBUTION NUMBERS AND THE ASSOCIATED FARM DEFINITIONS 
C 1 = UNIFORM ON THE INTERVAL RPARMl TO RPARM2 
C 2 = NORMAL WITH MEAN RPARMl AND VARIANCE RPARM2 
C 3 = EXPONENTIAL WITH RATE RPARMl 
C 4 = COUCHY WITH A = RPARMl AND B = RPARM2 
C 5 = GAMMA WITH SHAPE RPARM2 AND RATE RPARMl 
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C 6 = POISSON WITH RATE RPARMl 
C 7 = GEOMETRIC WITH P = RPARMl 
C 

REAL RPARMl, RPARM2,X 
INTEGER DISTN,SEED,IPARM,N 
REAL TEMP,VARIAT(1) 

IF(DISTN.LE.O .OR. DISTN.GE.8)THEN 
WRITE(10,*)’DISTN DOES NOT EXIT’ 

STOP 

ENDIF 

♦ 

GOTO (10,20,30,40,50,60,70),DISTN 

K< 

♦ GENERATE A UNIFORM BETWEEN RPARMl AND RPARM2 

* 

10 CONTINUE 

IF(RPARM1-RPARM2.EQ.0)THEN 
WRITE(10,*) ’ILLEGAL EQUAL RPARMS IN RANNUM’ 

STOP 

ENDIF 

IF(RPARM1 . GT.RPARM2)THEN 
TEMP = RPARMl 
RPARMl = RPARM2 
RPARM2 = TEMP 
ENDIF 

CALL LRND(SEED,VARIAT, 1,1,0) 

VARIAT(l) = RPARMl + (RPARM2-RPARM1)*VARIAT(1) 

GOTO 99 

* 

* GENERATE A NORMAL WITH MEAN RPARMl AND STD. DEV. RPARM2 

♦ 

20 CALL LNORM(SEED,VARIAT, 1,1,0) 

WRITE(10,’")’NORMAL(0,1) ’,VARIAT(1) 

VARIAT(l) = (VARIAT(1)’^RPARM2) + RPARMl 
GOTO 99 

* 

* GENERATE AN EXPONENTIAL WITH RATE RPARMl (1/MEAN) 

30 CONTINUE 

IF(RPARM1 . EQ.0.0)THEN 
WRITE(10,*)’ILLEGAL ZERO RATE IN RANNUM’ 

STOP 

ENDIF 
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CALL LEXPN(SEED, VARIAT, 1,1,0) 

VARIAT(l) = VARIAT(1)/RPARM1 
GOTO 99 

GENERATE A COUCHY WITH A = RPARMl AND B = RPARM2 

40 CONTINUE 

IF(RPARM2 .LE.0.0)THEN 

WRITE(10,*)’ILLEGAL COUCHY SPREAD IN RANNUM, B= ’,RPARM2 
STOP 
ENDIF 

CALL LCCHY(SEED, VARIAT, 1,1,0) 

VARIAT(l) = (VARIAT(1)*RPARM2) + RPARMl 
GOTO 99 

GENERATE A GAMMA WITH SHAPE RPARM2 AND RATE RPARMl 

50 CONTINUE 

IF(RPARM1 .LE.0.0)THEN 

WRITE(10,*) ’ILLEGAL NONPOSmVE GAMMA RATE IN RANNUM’ 
STOP 
ENDIF 

IF(RPARM2.LE.0.0)THEN 

WRITE(10,*)’ILLEGAL SHAPE PARAMETER IN RANNUM’ 

STOP 

ENDIF 

CALL LGAMA(SEED, VARIAT, 1,1, 0,RPARM2) 

VARIAT(l) = VARIAT(1)*(1.0/RPARM1) 

GOTO 99 

GENERATE A POISSON WITH RATE RPARMl 

60 CONTINUE 

IF(RPARM1 .LE.0.0)THEN 

WRITE(10,*)’ILLEGAL POISSON RATE IN RANNUM’ 

STOP 

ENDIF 

CALL LPOIS(SEED, VARIAT, 1,1,0,RPARM1) 

GOTO 99 

GENERATE A GEOMETRIC WITH P = RPARMl 



70 CONTINUE 



IF(RPAEM1 .LE.0.0)1HEN 

WRITE(10,*)’ILLEGAL GEOM. PROB. IN RANNUM’ 
STOP 
ENDIF 

CALL LGEOM(SEED,VARIAT,1,1,0,RPARM1) 

GOTO 99 

99 CONTINUE 
X = VARIAT(l) 

RETURN 

END 
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PROGRAM MOMENT 



C 

C TEAS PROGRAM USES TO COMPUTE TEIE SECOND AND THIRD 
C STANDARDIZED MOMENTS THAT SUPPORTS MOMENT MATCHING IN 

C GAMS PROGRAM. 

♦ 

* INITIALIZE DATA 

♦ 

INTEGER LIMIT, N,N1, 1, J, COUNT 

REAL UF, NEW(2), REPAGE, A, B, C, EX, EXSQ, EXC, EX3 
REAL STD, SECOND, THIRD, Cl, C2, AAC, OBJ 
REAL RAWDAT(1000,2), FB(IOOO), P(IOOO), TOT, TCOST 
DATA C1,C2,TOT,TCOST,AAC/100.0,500.0,0.0,0.0,0.0/ 

OPEN(UNIT = 1,FILE =’ LIFETIME’, STATUS =’OLD’) 

OPEN(UNIT = 2,FILE =’OPTAGE’, STATUS =’OLD’) 

OPEN(UNIT = 3,FILE =’RAWDATA’, STATUS =’OLD’) 

OPEN(UNIT = 4,FILE = ’MOMENT’) 

OPEN(UNIT = 5,FILE =’ COUNT’, STATUS =’OLD’) 

READ(1,*)UF 
READ(2,’")REPAGE 
READ(5,’^)COUNT 
PRINT’", ’COUNT= ’, COUNT 
N1 = 25 + COUNT 
DO 5 I = 1,1000 
P(I) = 0.0 
FB(I) = 0.0 
5 CONTINUE 
DO 44 I = 1,N1-1 

READ(3, 1 1)RAWDAT(1, 1) ,RAWDATa,2) 

11 F0RMAT(1X,F7.4,4X,F3.1) 

44 CONTINUE 
CLOSE(3) 

OPEN(UNIT=3,FILE= ’RAWDATA’) 

* 

* CALCULATION 

♦ 

IF(LIF.LE.REPAGE)THEN 
NEW(l) = UF 
NEW(2) = 1.0 
ELSE 

NEW(l) = REPAGE 
NEW(2) = 0.0 
ENDIF 
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CALL INSERT(N1,NEW,RAWDAT) 

DO 55 I = 1,N1 

WRITE(3,33) RAWDAT(I,1),RAWDAT(I,2) 

33 F0RMAT(1X,F7.4,4X,F3.1) 

55 CONTINUE 

IF(RAWDAT(1,2) .NE. 0.0)THEN 
FB(1) = (N1-1.0)/N1 
ELSE 

FB(1) = 1.0 
ENDIF 

DO 10 I = 2,N1 
IF(RAWDAT(I,2).EQ.0.0)THEN 
FB(I) = FB(I-1) 

ELSEIF(I.LT.N1)THEN 

FB(I) = ((N1-I)/(N1 + 1.0-I))*FB(I-1) 

ELSE 

FB(N1) = 0.0 
ENDIF 

10 CONTINUE 

P(l) = l.O-FB(l) 

DO 22 I = 2,N1 
P(I) = FB(I-1)-FB(I) 

22 CONTINUE 
EX = 0.0 
EXSQ = 0.0 
EXC = 0.0 
EX3 = 0.0 
DO 15 J = 1,N1 
TOT = TOT + RAWDAT(J,1) 

TCOST = TCOST + C2*RAWDAT(J,2) + C1*(1-RAWDAT(J,2)) 
A = RAWDAT(J,1)*P(J) 

B = (RAWDAT(J,1)**2)*P(J) 

C = (RAWDAT(J,1)**3)*P(J) 

EX = EX + A 
EXSQ = EXSQ + B 
EXC = EXC + C 
15 CONTINUE 

AAC = TCOST/TOT 

STD = (EXSQ - EX**2)**0.5 

EX3 = EXC + (2*EX**3) - (3*EX*EXSQ) 

SECOND = STD/EX 
THIRD = EX3/(STD**3) 
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* OUTPUT 

♦ 

WRrrE(4, 100)SECOND, THIRD, AAC, EX 
100 FORMAT(2X,F11.8,4X,F11.8,4X,F8.3,4X,F6.3) 

* CLOSE(5) 

* OPEN(UNIT = 5, FILE =’ COUNT’, STATUS =’OLD’) 

* WRITE(5,*)COUNT+l 
STOP 

END 

SUBROUTINE INSERT(UMIT,NEW,X) 

C THIS SUBROUTINE USES TO INSERT A PAIR ’NEW’ DATA INTO THE 
C DATA FILE ’X’ IN ASCENDING ORDER. 

INTEGER UMIT,J,K 
REAL NEW(2),X(1000,2) 

LOGICAL DONE 
DONE = .FALSE. 

J = 1 

5 IF((J.LT.LIMn).AND.(.NOT.DONE))THEN 
IF(X(J, 1).LT.NEW(1))THEN 
J =J+1 
ELSE 

DONE = .TRUE. 

ENDIF 
GOTO 5 
ENDIF 

IF(J. EQ.LIMIT)THEN 
X(J,1) = NEW(l) 

X(J,2) = NEW(2) 

ELSE 

IF(J.LT.LIMIT)THEN 
DO lOK = UMIT,J+1,-1 
X(K,1) = X(K-1,1) 

X(K,2) = X(K-1,2) 

10 CONTINUE 

X(J,1) = NEW(l) 

X(J,2) = NEW(2) 

ENDIF 

ENDIF 

RETURN 

END 
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PROGRAM DATA 



C 

C TfflS PROGRAM PROVTOES THE RESULT SEQUENTIAL ESTIMATED 
C ENTRIES TRANSITION RATES T(I,J) AND nSTTIAL STATE PROB. 

C ALPAH(I),AND ALSO ADJUST, PREPARE ALL THE INITIAL DATA THAT 
C WILL BE USED FOR THE NEXT SEQUENTIAL ESTIMATION. 

C THE VARIABLE ARE AS FOLLOWS: 

C A,B,C 

C D,E,F = THE ESTIMATED TRANSITION MATRIX 
C G,H,I 

C AL1,AL2,AL3 = THE ESTIMATED INITIAL STATES PROBABILITIES 
C X = RANDOM VARIABLE OF LIFETIME 
C OPTAGE = OPTIMAL AGE REPLACEMENT 
C SECOND = THE SECOND STANDARDIZED MOMENT 
C THIRD = THE THIRD STANDARDIZED MOMENT 
C OBJ = OBJECTIVE FUNCTION VALUE OF GAMS PROGRAM 

iK 

* INPUT 

INTEGER COUNT 

REAL A,B,C,D,E,F,G,H,I,AL1,AL2,AL3,X 

REAL EX, OPTCS, OPTAGE, Z, SECOND, THIRD, DI, OBJ, AAC, PI 

OPEN(UNIT = 10, FILE = ’OBJ’, STATUS =’OLD’) 

OPEN(UNIT = 11, FILE =’ ALPHA’, STATUS =’OLD’) 

OPEN(UNIT = 12, FILE = ’BRAVO’, STATUS =’OLD’) 

OPEN(UNIT = 13, FILE =’CHAREE’, STATUS =’OLD’) 

OPEN(UNIT = 14,FILE =’ DELTA’, STATUS =’OLD’) 

OPEN(UNIT = 15,FILE =’ ECHO’, STATUS =’OLD’) 

OPEN(UNIT = 16,FILE = ’FOXTROT’, STATUS =’OLD’) 

OPEN(UNIT = 17, FILE =’ GOLF’, STATUS =’OLD’) 

OPEN(UNIT = 18, FILE = ’HOTEL’, STATUS =’OLD’) 

OPEN(UNIT = 19,FILE = ’INDIA’, STATUS =’OLD’) 

OPEN(UNIT = 20,FILE =’AL1’,STATUS =’OLD’) 

OPEN(UNIT = 21, FILE =’AL2’, STATUS =’OLD’) 

OPEN(UNIT = 22, FILE =’AL3’,STATUS =’OLD’) 

OPEN(UNIT = 23, FILE =’ OPTAGE’, STATUS =’OLD’) 

OPEN(UNIT = 24,FILE = ’LIFETIME’, STATUS =’OLD’) 

OPEN(UNIT = 25,FILE = ’MOMENT’, STATUS =’OLD’) 

OPEN(UNIT = 26,FILE =’ COUNT’, STATUS =’OLD’) 

OPEN(UNIT = 27, FILE =’OPTCOST’, STATUS =’OLD’) 

CALL EXCMS(’FILEDEF 50 DISK FILE ESTDATl (DISP MOD’) 

CALL EXCMS(’FILEDEF 60 DISK FILE ESTDAT2 (DISP MOD’) 

CALL EXCMS(’FILEDEF 70 DISK FILE ESTDAT3 (DISP MOD’) 
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OPEN(UNIT = 30, FILE = ’SCALAR!’) 
READ(10,’^)OBJ 
IF(OBJ .GT. 0.05)THEN 
A = -1.0 
B = 0.3 
C = 0.3 
D = 0.3 
E = -1.0 
F = 0.3 
G = 0.3 
H = 0.3 
I = -1.0 
ELSE 

READ(11,*)A 

READ(12,*)B 

READ(13,*)C 

READ(14,*)D 

READ(15,*)E 

READ(16,*)F 

READ(17,*)G 

READ(18,*)H 

READ(19,*)I 

ENDIF 

READ(20,*)AL1 
READ(21,*)AL2 
READ(22,’^)AL3 
READ(23, ’^)OPTAGE,PI 
READ(24,’^)X 

READ(25, ’^)SECOND, THIRD, AAC, EX 

READ(26,’^)COUNT 

READ(27,’^)OPTCS 

CLOSE(26) 

PRINT*, ’OBJ=’,OBJ 



IF(OPTAGE .GT. X)THEN 
Z = X 
DI = 1.0 
ELSE 

Z = OPTAGE 
DI = 0.0 
ENDIF 
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* OUTPUT 
» 

OPEN(UNIT = 26, FILE =’ COUNT’, STATUS =’OLD’) 
WRITE(26, *)COUNT+ 1 
WRITE(50,200)COUNT,A,B,C,D,E,F,G,H,I 
WRITE(60,300)COUNT,OPTCS,OPTAGE,X,Z,DI,AAC 
WRITE(70,100)COUNT,AL1,AL2,AL3,SECOND,THIRD,OBJ,PI 



PREPARE THE DATA FOR "GAMS" PROGRAM 



WRITE(30,*)’ 
WRITE(30,400)’ 
WRITE(30,400)’ 
WRITE(30,400)’ 
WRITE(30,400)’ 
WRITE(30,400)’ 
WRITE(30,400)’ 
WRITE(30,400)’ 
WRITE(30,400)’ 
WRITE(30,400)’ 
WRITE(30,500)’ 
WRITE(30,500)’ 
WRITE(30,500)’ 
WRITE(30,600)’ 
WRITE(30,600)’ 
WRITE(30,700)’ 



SCALAR’ 

AA INITIAL SOLUTION OF A/’, A,’/’ 

BB INITIAL SOLUTION OF B /’,B,’ /’ 

CC INITIAL SOLUTION OF C /’,C,’ /’ 

DD INITIAL SOLUTION OF D /’,D,’ /’ 

EE INITIAL SOLUTION OF A/’, E,’ /’ 

FF INITIAL SOLUTION OF B/’,F,’/’ 

GG INITIAL SOLUTION OF C/’,G,’/’ 

HH INITIAL SOLUTION OF D /’,H,’ /’ 
n INITIAL SOLUTION OF D /’ 

ALPl INITIAL SOLUTION OF ALl /’,AL1,’ /’ 

ALP2 INITIAL SOLUTION OF AL2 /’,AL2,’ /’ 

ALP3 INITIAL SOLUTION OF AL3 /’,AL3,’ /’ 

MM2 SECOND STANDARD MOMENT /’, SECOND,’/’ 
MM3 THIRD STANDARD MOMENT /’,THIRD,’/’ 
EX EXPECTED VALUE OF LIFETIME /’,EX,’ / ;’ 
100 FORMAT(lX,I3,3(2X,F5.2),2(2X,F8.5),2X,F4.2,2X,E11.4) 

200 F0RMAT(1X,I3,2X,9(2X,F6.2)) 

300 FORMAT(lX,I3,2X,F8.3,2X,F5.2,2(4X,F7.4),4X,F3.1,4X,F7.3) 

400 FORMAT(A35,F6.2,A2) 

500 FORMAT(A39,F5.2,A2) 

600 F0RMAT(A36,F11.8,A1) 

700 FORMAT(A40,F6.3,A4) 

STOP 

END 
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